################
#Public Opinion#
################
#This to obtain:
#table_4.tex
#table_A13.tex
#table_A15.tex

#Figure_A17a.jpg
#Figure_A17b.jpg
#Figure_A18a.jpg
#Figure_A18b.jpg
#Figure_A18c.jpg
#Figure_A18d.jpg
#Figure_A18e.jpg
#Figure_A18f.jpg
#Figure_A18g.jpg
#Figure_A18h.jpg
#Figure_A19a.jpg
#Figure_A19b.jpg
#Figure_A19c.jpg
#Figure_A19d.jpg
#Figure_A19e.jpg
#Figure_A19f.jpg
#Figure_A19g.jpg
#Figure_A19h.jpg



rm(list = ls())

library("stargazer")
library("sandwich")
library("Greg")
library("broom")
library("ggplot2")
library("stringr")
library("multiwayvcov")
library("stringi")
library("readxl")
library("sf")
library("lemon")

#Loading up the dataset
setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#Functions for coefficient plots

rerun_regression<-function(regression_name, string_dv){
  #Process:
  #-takes already existing regression objects, associates them with an 
  #"official" label, calculates robust SE and CI and puts them in a tibble
  
  #Input:
  #-Regression model object after lm
  #-The label for the DV that you want displayed in ggplot
  
  #Output:
  #-A tibble with regression coefficients from the model that you ran with:
  #--robust SE
  #--robust CI
  
  #Recreate dataframe from the regression object
  remake_data<-regression_name$model
  #Scale the DV
  remake_data$dv_scaled<-scale(remake_data[1])
  #Take the string of the DV and IV
  dv<-regression_name$terms[[2]]
  ivs<-regression_name$terms[[3]]
  ivs_str<-paste(ivs)[2]
  dv_str<-paste(dv, "~", collapse = "")
  dv_ivs_str<-paste(dv_str, ivs_str, sep=" ")
  dv_str_scaled<-paste("dv_scaled", "~", collapse = "")
  dvscales_ivs_str<-paste(dv_str_scaled, ivs_str, sep=" ")
  #Create specification and run model
  specification = as.formula(dvscales_ivs_str)
  model <- lm(specification, data = remake_data)
  no_obs<-nobs(model)
  no_obs<-paste("Obs: ", as.character(no_obs), sep = "")
  #Calculate robust SE
  model_cov <- vcovHC(model, type = "HC")
  model_se <- sqrt(diag(model_cov))
  #Create tibble
  model_tidy<-tidy(model)
  model_tidy$no_obs<-no_obs
  
  model_se_tidy<-tidy(model_se)
  names(model_se_tidy)[1] <- "term"
  names(model_se_tidy)[2] <- "robust.se"
  model_tidy2<-left_join(model_tidy, model_se_tidy, by = c("term"="term"))
  
  #Calculate CI
  df_ci<-as.data.frame(confint_robust(model, level = 0.95, HC_type = "HC", t_distribution = FALSE))
  df_ci <- cbind(term = rownames(df_ci), df_ci)
  names(df_ci)[2] <- "lb_95"
  names(df_ci)[3] <- "ub_95"
  
  #Merge CI to the dataframe
  model_tidy3<-left_join(model_tidy2, df_ci, by = c("term"="term"))
  
  #Only keep the treat variable
  model_tidy4<-subset(model_tidy3, term=="treat")
  model_tidy4$dv_official<-string_dv
  
  return(model_tidy4)
}



return_graphs_coef<-function(models_run, models_run_alias){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  #return(big_data)
  mean_feature <- ggplot(data = big_data, aes(x=estimate, y=as.numeric(reorder(dv_official, desc(counter))))) +
    geom_point(size = 2) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs1),
                       labels = irislabs1,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    scale_x_continuous(name = "Estimate") +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=16))  
  #return(big_data)
  return(mean_feature)
}


return_graphs_coef_urb_rur<-function(models_run, models_run_alias){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  
  big_data$model2<-NA
  big_data$model2[stri_detect_fixed(big_data$dv_official, "Urban")]<-"Urban"
  big_data$model2[stri_detect_fixed(big_data$dv_official, "Rural")]<-"Rural"
  big_data$model2[is.na(big_data$model2)]<-"All"
  
  #define as factor in the desired order
  big_data$model2<-factor(big_data$model2, levels=c("All", "Urban", "Rural"))
  
  big_data$model<-big_data$dv_official
  big_data$model[stri_detect_fixed(big_data$dv_official, "Urban")]<-""
  big_data$model[stri_detect_fixed(big_data$dv_official, "Rural")]<-""
  
  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  irislabs3 <- rev(big_data$model)
  
  
  mean_feature<-ggplot(data = big_data, aes(x=estimate, 
                                            y=as.numeric(reorder(dv_official, desc(counter))),
                                            colour= model2))+
    geom_point(aes(shape=model2), position = position_dodge(0.5), size = 4) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0,
                  position=position_dodge(width=0.5))+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs3),
                       labels = irislabs3,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    scale_colour_manual(name = "Sample", 
                        label = c("All", "Urban", "Rural"),
                        values=c("grey5", "grey50", "grey80"))+
    scale_shape_manual(name = "Sample", 
                       label = c("All", "Urban", "Rural"),
                       values=c(16, 15, 17))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    #scale_y_discrete(name = "") +
    scale_x_continuous(name = "Estimate") +
    #ggtitle(glue("{title_spec}")) +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=16),
          legend.title=element_text(size=14), 
          legend.text=element_text(size=14))
  #return(big_data)
  return(mean_feature)
}


return_graphs_coef_cohort<-function(models_run, models_run_alias){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  
  big_data$model2<-NA
  big_data$model2[stri_detect_fixed(big_data$dv_official, "1945")]<-"Born Before 1945"
  big_data$model2[stri_detect_fixed(big_data$dv_official, "1945-1990")]<-"Born 1945-1990"
  big_data$model2[stri_detect_fixed(big_data$dv_official, "1990-2016")]<-"Born 1990-2016"
  
  
  #define as factor in the desired order
  big_data$model2<-factor(big_data$model2, levels=c("Born Before 1945", "Born 1945-1990", "Born 1990-2016"))
  
  big_data$model<-big_data$dv_official
  big_data$model[stri_detect_fixed(big_data$dv_official, "1945-1990")]<-""
  big_data$model[stri_detect_fixed(big_data$dv_official, "1990-2016")]<-""
  big_data$model<-str_replace(big_data$model, "Cohort Before 1945\n ", "")
  

  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  irislabs3 <- rev(big_data$model)
  
  
  #return(big_data)
  mean_feature<-ggplot(data = big_data, aes(x=estimate, 
                                            y=as.numeric(reorder(dv_official, desc(counter))),
                                            colour= model2))+
    geom_point(aes(shape=model2), position = position_dodge(0.5), size = 4) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0,
                  position=position_dodge(width=0.5))+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs3),
                       labels = irislabs3,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    scale_colour_manual(name = "Sample", 
                        label = c("Born Before 1945", "Born 1945-1990", "Born 1990-2016"),
                        values=c("grey5", "grey50", "grey80"))+
    scale_shape_manual(name = "Sample", 
                       label = c("Born Before 1945", "Born 1945-1990", "Born 1990-2016"),
                       values=c(16, 15, 17))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    scale_x_continuous(name = "Estimate") +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=14),
          axis.text.y = element_text(size=14),
          axis.title=element_text(size=16),
          legend.title=element_text(size=14), 
          legend.text=element_text(size=14))
  #return(big_data)
  return(mean_feature)
}


clean_col_labels<-function(column_labels){
  #Process:
  #-takes list of column labels for stargazer
  #and removes the unnecessary tags
  
  #Input:
  #-list of column labels for stargazer
  
  #Output:
  #-retruns list of clean labels
  df<-list()
  length(df)<-length(column_labels)
  counter = 1
  for (i in column_labels){
    #print(i)
    j<-str_remove(i, '\\{')
    k<-str_remove(j, '\\}')
    l<-str_remove(k, '\\\\')
    m<-str_remove(l, 'shortstack')
    n<-gsub("\\\\", "\n", m)
    o<-gsub("\n\n", "\n", n)
    p<-gsub(",", "", o)
    df[counter] = p
    counter = counter +1
    df<-unlist(df)
  }
  return(df)
}

lits_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/data_lits.xlsx", 
                       sheet=1, col_names = TRUE, skip = 0, guess_max = 10000)

data_before_1945<-subset(lits_data, year_birth<1945)
data_before_1945_1990<-subset(lits_data, year_birth<1990 & year_birth>1945)
data_before_1990_2016<-subset(lits_data, year_birth>1990)
data_only_2010<-subset(lits_data, survey_2010==1)
data_only_2016<-subset(lits_data, survey_2016==1)
data20km<-subset(lits_data, krajna6_distance<20)
data40km<-subset(lits_data, krajna6_distance<40)
data60km<-subset(lits_data, krajna6_distance<60)


#MAPS
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")
HRV_adm0<-st_simplify(HRV_adm0,  dTolerance = 0.005)

mil_col_short_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")
mil_col_short_polygon<-st_simplify(mil_col_short_polygon,  dTolerance = 0.005)

#Indentifying unique coordinates
data <- st_as_sf(x = lits_data,                         
                  coords = c("longitude", "latitude"),
                  crs = st_crs(HRV_adm0))

data$x_lon<-sf::st_coordinates(st_centroid(data))[,1]
data$y_lat<-sf::st_coordinates(st_centroid(data))[,2]
data2<-subset(data, !duplicated(data[c("y_lat", "x_lon")]))

########
#MAP 1 #
########

col = c("2006"="grey40",
        "2010"="grey40",
        "2016" = "grey40")

order_col<-c("2006", "2010", "2016")

lits_map<-ggplot()+
  geom_sf(data = HRV_adm0, 
          colour = "black",
          fill = NA, linewidth = 0.8)+
  geom_point(data=data2, 
             size = 2,
             aes(x = x_lon, y= y_lat, 
                 color=factor(survey_year),
                 shape = factor(survey_year)), 
             alpha=0.5,
             position=position_jitter(width=.05, height=.05))+
  scale_colour_manual(values = col,
                      name = "Survey Year",
                      breaks = order_col)+
  scale_shape_manual(values = c(3, 16, 17))+
  guides(colour = guide_legend("Survey Year"), 
         size = guide_legend("Survey Year"),
         shape = guide_legend("Survey Year"))+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 1.5)+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(legend.position="left")+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                 #Above you mention how long the scalebar should be
                 y.min = 46.4, y.max = 46.7,
                 #Above you mention how thick the scalebar should be
                 dist = 50, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.2,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)

lits_map2<-reposition_legend(lits_map, 'bottom left')

ggsave(lits_map2,file="./Dropbox/Legacies_Project/Paper/figures/Figure_A17a.jpg", 
       height=15, width=15, units = "cm", dpi=300)


data2$dist_group<-NA
data2$dist_group[data2$krajna6_distance<=20]<-"20km"
data2$dist_group[data2$krajna6_distance>20 & data2$krajna6_distance<=40]<-"40km"
data2$dist_group[data2$krajna6_distance>40 & data2$krajna6_distance<=120]<-"60km"

col = c("20km"="red",
        "40km"="yellow",
        "60km" = "green")

order_col<-c("20km", "40km", "60km")

lits_map<-ggplot()+
  geom_sf(data = HRV_adm0, 
          colour = "black",
          fill = NA, lindwidth = 0.8)+
  geom_point(data = data2, aes(x = x_lon, y= y_lat,
                               color = dist_group),
             alpha=0.5,
             position=position_jitter(width=.08, height=.08))+
  scale_colour_manual(values=col, 
                      guide = guide_legend(title = "Distance Border",
                                           override.aes = list(size=3)))+
  geom_sf(data = mil_col_short_polygon, fill = NA, colour = "blue", lindwidth = 1.5)+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(legend.position="left")+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                 #Above you mention how long the scalebar should be
                 y.min = 46.4, y.max = 46.7,
                 #Above you mention how thick the scalebar should be
                 dist = 50, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.2,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)

lits_map2<-reposition_legend(lits_map, 'bottom left')

ggsave(lits_map2,file="./Dropbox/Legacies_Project/Paper/figures/Figure_A17b.jpg", 
       height=15, width=15, units = "cm", dpi=300)




#Models
trust_people<-lm(trust_people~treat+
                POINT_X+POINT_Y+zagreb_distance+
                bfe1+bfe2+bfe3+
                bfe4+bfe5+bfe6+
                bfe7+
                age_group +income +edu+
                survey_2006+
                survey_2010+
                survey_2016,
              data=data, na.action=na.exclude)
summary(trust_people)
trust_people_cov <- vcovHC(trust_people, type = "HC")
trust_people_se <- sqrt(diag(trust_people_cov))


trust_people_20km<-lm(trust_people~treat+
                            POINT_X+POINT_Y+zagreb_distance+
                            bfe1+bfe2+bfe3+
                            bfe4+bfe5+bfe6+
                            bfe7+
                            age_group +income +edu+
                            survey_2006+
                            survey_2010+
                            survey_2016,
                          data=subset(data, krajna6_distance<20),
                          na.action=na.exclude)
summary(trust_people_20km)


trust_people_40km<-lm(trust_people~treat+
                            POINT_X+POINT_Y+zagreb_distance+
                            bfe1+bfe2+bfe3+
                            bfe4+bfe5+bfe6+
                            bfe7+
                            age_group +income +edu+
                            survey_2006+
                            survey_2010+
                            survey_2016,
                          data=subset(data, krajna6_distance<40),
                          na.action=na.exclude)
summary(trust_people_40km)


trust_people_60km<-lm(trust_people~treat+
                            POINT_X+POINT_Y+zagreb_distance+
                            bfe1+bfe2+bfe3+
                            bfe4+bfe5+bfe6+
                            bfe7+
                            age_group +income +edu+
                            survey_2006+
                            survey_2010+
                            survey_2016,
                          data=subset(data, krajna6_distance<60),
                          na.action=na.exclude)
summary(trust_people_60km)


trust_people_b1945<-lm(scale(trust_people)~treat+
                               POINT_X+POINT_Y+zagreb_distance+
                               bfe1+bfe2+bfe3+
                               bfe4+bfe5+bfe6+
                               bfe7+
                               age_group +income +edu+
                               survey_2006+
                               survey_2010+
                               survey_2016,
                             data=data_before_1945, na.action=na.exclude)
summary(trust_people_b1945)
trust_people_b1945_cov <- vcovHC(trust_people_b1945, type = "HC")
trust_people_b1945_se <- sqrt(diag(trust_people_b1945_cov))


trust_people_1945_1990<-lm(scale(trust_people)~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=data_before_1945_1990, na.action=na.exclude)
summary(trust_people_1945_1990)
trust_people_1945_1990_cov <- vcovHC(trust_people_1945_1990, type = "HC")
trust_people_1945_1990_se <- sqrt(diag(trust_people_1945_1990_cov))


trust_people_1990_2016<-lm(scale(trust_people)~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=data_before_1990_2016, na.action=na.exclude)
summary(trust_people_1990_2016)
trust_people_1990_2016_cov <- vcovHC(trust_people_1945_1990, type = "HC")
trust_people_1990_2016_se <- sqrt(diag(trust_people_1990_2016_cov))


##################
#Trust Presidency#
##################

trust_presidency_20km<-lm(trust_presidency~treat+
                       POINT_X+POINT_Y+zagreb_distance+
                       bfe1+bfe2+bfe3+
                       bfe4+bfe5+bfe6+
                       bfe7+
                       age_group +income +edu+
                       survey_2006+
                       survey_2010+
                       survey_2016,
                     data=subset(data, krajna6_distance<20),
                     na.action=na.exclude)
summary(trust_presidency_20km)


trust_presidency_40km<-lm(trust_presidency~treat+
                            POINT_X+POINT_Y+zagreb_distance+
                            bfe1+bfe2+bfe3+
                            bfe4+bfe5+bfe6+
                            bfe7+
                            age_group +income +edu+
                            survey_2006+
                            survey_2010+
                            survey_2016,
                          data=subset(data, krajna6_distance<40),
                          na.action=na.exclude)
summary(trust_presidency_40km)


trust_presidency_60km<-lm(trust_presidency~treat+
                            POINT_X+POINT_Y+zagreb_distance+
                            bfe1+bfe2+bfe3+
                            bfe4+bfe5+bfe6+
                            bfe7+
                            age_group +income +edu+
                            survey_2006+
                            survey_2010+
                            survey_2016,
                          data=subset(data, krajna6_distance<60),
                          na.action=na.exclude)
summary(trust_presidency_60km)


trust_presidency<-lm(trust_presidency~treat+
                       POINT_X+POINT_Y+zagreb_distance+
                       bfe1+bfe2+bfe3+
                       bfe4+bfe5+bfe6+
                       bfe7+
                       age_group +income +edu+
                       survey_2006+
                       survey_2010+
                       survey_2016,
                     data=subset(data),
                     na.action=na.exclude)
summary(trust_presidency)
trust_presidency_cov <- vcovHC(trust_presidency, type = "HC")
trust_presidency_se <- sqrt(diag(trust_presidency_cov))


trust_presidency_b1945<-lm(trust_presidency~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=data_before_1945,
                           na.action=na.exclude)
summary(trust_presidency_b1945)
trust_presidency_b1945_cov <- vcovHC(trust_presidency_b1945, type = "HC")
trust_presidency_b1945_se <- sqrt(diag(trust_presidency_b1945_cov))

trust_presidency_1945_1990<-lm(trust_presidency~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=data_before_1945_1990,
                           na.action=na.exclude)
summary(trust_presidency_1945_1990)
trust_presidency_1945_1990_cov <- vcovHC(trust_presidency_1945_1990, type = "HC")
trust_presidency_1945_1990_se <- sqrt(diag(trust_presidency_1945_1990_cov))


trust_presidency_1990_2016<-lm(trust_presidency~treat+
                                 POINT_X+POINT_Y+zagreb_distance+
                                 bfe1+bfe2+bfe3+
                                 bfe4+bfe5+bfe6+
                                 bfe7+
                                 age_group +income +edu+
                                 survey_2006+
                                 survey_2010+
                                 survey_2016,
                               data=data_before_1990_2016,
                               na.action=na.exclude)
summary(trust_presidency_1990_2016)
trust_presidency_1990_2016_cov <- vcovHC(trust_presidency_1990_2016, type = "HC")
trust_presidency_1990_2016_se <- sqrt(diag(trust_presidency_1990_2016_cov))

stargazer(trust_presidency_b1945, 
          trust_presidency_1945_1990,
          trust_presidency_1990_2016,
          type= "text", keep = c("treat"),
          se = list(trust_presidency_b1945_se,
                    trust_presidency_1945_1990_se,
                    trust_presidency_1990_2016_se))




###############
# trust_family#
###############

trust_family_b1945<-lm(trust_family~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=data_before_1945, na.action = na.exclude)
summary(trust_family_b1945)
trust_family_b1945_cov <- vcovHC(trust_family_b1945, type = "HC")
trust_family_b1945_se <- sqrt(diag(trust_family_b1945_cov))

trust_family_1945_1990<-lm(trust_family~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=data_before_1945_1990, na.action = na.exclude)
summary(trust_family_1945_1990)
trust_family_1945_1990_cov <- vcovHC(trust_family_1945_1990, type = "HC")
trust_family_1945_1990_se <- sqrt(diag(trust_family_1945_1990_cov))


trust_family_1990_2016<-lm(trust_family~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=data_before_1990_2016, na.action = na.exclude)
summary(trust_family_1990_2016)
trust_family_1990_2016_cov <- vcovHC(trust_family_1990_2016, type = "HC")
trust_family_1990_2016_se <- sqrt(diag(trust_family_1990_2016_cov))


stargazer(trust_family_b1945, 
          trust_family_1945_1990,
          trust_family_1990_2016,
          type= "text", keep = c("treat"),
          se = list(trust_family_b1945_se,
                    trust_family_1945_1990_se,
                    trust_family_1990_2016_se))


trust_family<-lm(trust_family~treat+
                        POINT_X+POINT_Y+zagreb_distance+
                        bfe1+bfe2+bfe3+
                        bfe4+bfe5+bfe6+
                        bfe7+
                        age_group +income +edu+
                        survey_2006+
                        survey_2010+
                        survey_2016,
                      data=subset(data),
                      na.action = na.exclude)
summary(trust_family)
trust_family_cov <- vcovHC(trust_family, type = "HC")
trust_family_se <- sqrt(diag(trust_family_cov))



trust_family_20km<-lm(scale(trust_family)~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=subset(data, krajna6_distance<20),
                            na.action = na.exclude)
summary(trust_family_20km)

trust_family_40km<-lm(scale(trust_family)~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=subset(data, krajna6_distance<40),
                           na.action = na.exclude)
summary(trust_family_40km)

trust_family_60km<-lm(scale(trust_family)~treat+
                             POINT_X+POINT_Y+zagreb_distance+
                             bfe1+bfe2+bfe3+
                             bfe4+bfe5+bfe6+
                             bfe7+
                             age_group +income +edu+
                             survey_2006+
                             survey_2010+
                             survey_2016,
                           data=subset(data, krajna6_distance<60),
                           na.action = na.exclude)
summary(trust_family_60km)



#############
# trust_foreign_investors
############

trust_foreign_investors_20km<-lm(trust_foreign_investors~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=subset(data, krajna6_distance<20), na.action=na.exclude)
summary(trust_foreign_investors_20km)


trust_foreign_investors_40km<-lm(trust_foreign_investors~treat+
                                   POINT_X+POINT_Y+zagreb_distance+
                                   bfe1+bfe2+bfe3+
                                   bfe4+bfe5+bfe6+
                                   bfe7+
                                   age_group +income +edu+
                                   survey_2006+
                                   survey_2010+
                                   survey_2016,
                                 data=subset(data, krajna6_distance<40), na.action=na.exclude)
summary(trust_foreign_investors_40km)

trust_foreign_investors_60km<-lm(trust_foreign_investors~treat+
                                   POINT_X+POINT_Y+zagreb_distance+
                                   bfe1+bfe2+bfe3+
                                   bfe4+bfe5+bfe6+
                                   bfe7+
                                   age_group +income +edu+
                                   survey_2006+
                                   survey_2010+
                                   survey_2016,
                                 data=subset(data, krajna6_distance<60), na.action=na.exclude)
summary(trust_foreign_investors_60km)



trust_foreign_investors<-lm(trust_foreign_investors~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=subset(data), na.action=na.exclude)
summary(trust_foreign_investors)
trust_foreign_investors_cov <- vcovHC(trust_foreign_investors, type = "HC")
trust_foreign_investors_se <- sqrt(diag(trust_foreign_investors_cov))


trust_foreign_investors_b1945<-lm(trust_foreign_investors~treat+
                                    POINT_X+POINT_Y+zagreb_distance+
                                    bfe1+bfe2+bfe3+
                                    bfe4+bfe5+bfe6+
                                    bfe7+
                                    age_group +income +edu+
                                    survey_2006+
                                    survey_2010+
                                    survey_2016,
                                  data=data_before_1945, na.action=na.exclude)
summary(trust_foreign_investors_b1945)
trust_foreign_investors_b1945_cov <- vcovHC(trust_foreign_investors_b1945, type = "HC")
trust_foreign_investors_b1945_se <- sqrt(diag(trust_foreign_investors_b1945_cov))


trust_foreign_investors_1945_1990<-lm(trust_foreign_investors~treat+
                                    POINT_X+POINT_Y+zagreb_distance+
                                    bfe1+bfe2+bfe3+
                                    bfe4+bfe5+bfe6+
                                    bfe7+
                                    age_group +income +edu+
                                    survey_2006+
                                    survey_2010+
                                    survey_2016,
                                  data=data_before_1945_1990, na.action=na.exclude)
summary(trust_foreign_investors_1945_1990)
trust_foreign_investors_1945_1990_cov <- vcovHC(trust_foreign_investors_1945_1990, type = "HC")
trust_foreign_investors_1945_1990_se <- sqrt(diag(trust_foreign_investors_1945_1990_cov))


trust_foreign_investors_1990_2016<-lm(trust_foreign_investors~treat+
                                        POINT_X+POINT_Y+zagreb_distance+
                                        bfe1+bfe2+bfe3+
                                        bfe4+bfe5+bfe6+
                                        bfe7+
                                        age_group +income +edu+
                                        survey_2006+
                                        survey_2010+
                                        survey_2016,
                                      data=data_before_1945_1990, na.action=na.exclude)
summary(trust_foreign_investors_1990_2016)
trust_foreign_investors_1990_2016_cov <- vcovHC(trust_foreign_investors_1990_2016, type = "HC")
trust_foreign_investors_1990_2016_se <- sqrt(diag(trust_foreign_investors_1990_2016_cov))


#####################
#bribery_road_police#
#####################

bribery_road_police_20km<-lm(bribery_road_police~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=subset(data, krajna6_distance<20), na.action = na.exclude)
summary(bribery_road_police_20km)

bribery_road_police_40km<-lm(bribery_road_police~treat+
                               POINT_X+POINT_Y+zagreb_distance+
                               bfe1+bfe2+bfe3+
                               bfe4+bfe5+bfe6+
                               bfe7+
                               age_group +income +edu+
                               survey_2006+
                               survey_2010+
                               survey_2016,
                             data=subset(data, krajna6_distance<40), na.action = na.exclude)
summary(bribery_road_police_40km)


bribery_road_police_60km<-lm(bribery_road_police~treat+
                               POINT_X+POINT_Y+zagreb_distance+
                               bfe1+bfe2+bfe3+
                               bfe4+bfe5+bfe6+
                               bfe7+
                               age_group +income +edu+
                               survey_2006+
                               survey_2010+
                               survey_2016,
                             data=subset(data, krajna6_distance<60), na.action = na.exclude)
summary(bribery_road_police_60km)

bribery_road_police<-lm(bribery_road_police~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                    data=subset(data), na.action = na.exclude)
summary(bribery_road_police)
bribery_road_police_cov <- vcovHC(bribery_road_police, type = "HC")
bribery_road_police_se <- sqrt(diag(bribery_road_police_cov))


bribery_road_police_b1945<-lm(bribery_road_police~treat+
                                POINT_X+POINT_Y+zagreb_distance+
                                bfe1+bfe2+bfe3+
                                bfe4+bfe5+bfe6+
                                bfe7+
                                age_group +income +edu+
                                survey_2006+
                                survey_2010+
                                survey_2016,
                              data=data_before_1945, na.action = na.exclude)
summary(bribery_road_police_b1945)
bribery_road_police_b1945_cov <- vcovHC(bribery_road_police_b1945, type = "HC")
bribery_road_police_b1945_se <- sqrt(diag(bribery_road_police_b1945_cov))


bribery_road_police_1945_1990<-lm(bribery_road_police~treat+
                                POINT_X+POINT_Y+zagreb_distance+
                                bfe1+bfe2+bfe3+
                                bfe4+bfe5+bfe6+
                                bfe7+
                                age_group +income +edu+
                                survey_2006+
                                survey_2010+
                                survey_2016,
                              data=data_before_1945_1990, na.action = na.exclude)
summary(bribery_road_police_1945_1990)
bribery_road_police_1945_1990_cov <- vcovHC(bribery_road_police_1945_1990, type = "HC")
bribery_road_police_1945_1990_se <- sqrt(diag(bribery_road_police_1945_1990_cov))


bribery_road_police_1990_2016<-lm(bribery_road_police~treat+
                                    POINT_X+POINT_Y+zagreb_distance+
                                    bfe1+bfe2+bfe3+
                                    bfe4+bfe5+bfe6+
                                    bfe7+
                                    age_group +income +edu+
                                    survey_2006+
                                    survey_2010+
                                    survey_2016,
                                  data=data_before_1990_2016, na.action = na.exclude)
summary(bribery_road_police_1990_2016)
bribery_road_police_1990_2016_cov <- vcovHC(bribery_road_police_1990_2016, type = "HC")
bribery_road_police_1990_2016_se <- sqrt(diag(bribery_road_police_1990_2016_cov))



################
#bribery_unempl#
################

bribery_unempl_20km<-lm(bribery_unempl~treat+
                     POINT_X+POINT_Y+zagreb_distance+
                     bfe1+bfe2+bfe3+
                     bfe4+bfe5+bfe6+
                     bfe7+
                     age_group +income +edu+
                     survey_2006+
                     survey_2010+
                     survey_2016,
                   data=subset(data, krajna6_distance<20),
                   na.action = na.exclude)
summary(bribery_unempl_20km)

bribery_unempl_40km<-lm(bribery_unempl~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=subset(data, krajna6_distance<40),
                        na.action = na.exclude)
summary(bribery_unempl_40km)

bribery_unempl_60km<-lm(bribery_unempl~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=subset(data, krajna6_distance<60),
                        na.action = na.exclude)
summary(bribery_unempl_60km)

bribery_unempl<-lm(bribery_unempl~treat+
                     POINT_X+POINT_Y+zagreb_distance+
                     bfe1+bfe2+bfe3+
                     bfe4+bfe5+bfe6+
                     bfe7+
                     age_group +income +edu+
                     survey_2006+
                     survey_2010+
                     survey_2016,
                   data=subset(data),
                   na.action = na.exclude)
summary(bribery_unempl)
bribery_unempl_cov <- vcovHC(bribery_unempl, type = "HC")
bribery_unempl_se <- sqrt(diag(bribery_unempl_cov))


bribery_unempl_b1945<-lm(bribery_unempl~treat+
                           POINT_X+POINT_Y+zagreb_distance+
                           bfe1+bfe2+bfe3+
                           bfe4+bfe5+bfe6+
                           bfe7+
                           age_group +income +edu+
                           survey_2006+
                           survey_2010+
                           survey_2016,
                         data=data_before_1945,
                         na.action = na.exclude)
summary(bribery_unempl_b1945)
bribery_unempl_b1945_cov <- vcovHC(bribery_unempl_b1945, type = "HC")
bribery_unempl_b1945_se <- sqrt(diag(bribery_unempl_b1945_cov))



bribery_unempl_1945_1990<-lm(bribery_unempl~treat+
                           POINT_X+POINT_Y+zagreb_distance+
                           bfe1+bfe2+bfe3+
                           bfe4+bfe5+bfe6+
                           bfe7+
                           age_group +income +edu+
                           survey_2006+
                           survey_2010+
                           survey_2016,
                         data=data_before_1945_1990,
                         na.action = na.exclude)
summary(bribery_unempl_1945_1990)
bribery_unempl_1945_1990_cov <- vcovHC(bribery_unempl_1945_1990, type = "HC")
bribery_unempl_1945_1990_se <- sqrt(diag(bribery_unempl_1945_1990_cov))



bribery_unempl_1990_2016<-lm(bribery_unempl~treat+
                               POINT_X+POINT_Y+zagreb_distance+
                               bfe1+bfe2+bfe3+
                               bfe4+bfe5+bfe6+
                               bfe7+
                               age_group +income +edu+
                               survey_2006+
                               survey_2010+
                               survey_2016,
                             data=data_before_1990_2016,
                             na.action = na.exclude)
summary(bribery_unempl_1990_2016)
bribery_unempl_1990_2016_cov <- vcovHC(bribery_unempl_1990_2016, type = "HC")
bribery_unempl_1990_2016_se <- sqrt(diag(bribery_unempl_1990_2016_cov))


stargazer(bribery_unempl_b1945, 
          bribery_unempl_1945_1990,
          bribery_unempl_1990_2016,
          keep = c("treat"),
          type = "text",
          covariate.labels=c("Military Territory"),
          se=list(bribery_unempl_b1945_se, 
                  bribery_unempl_1945_1990_se, 
                  bribery_unempl_1990_2016_se),
          dep.var.labels.include = F,
          omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE)


################
#bribery_documents#
################
bribery_documents<-lm(bribery_documents~treat+
                     POINT_X+POINT_Y+zagreb_distance+
                     bfe1+bfe2+bfe3+
                     bfe4+bfe5+bfe6+
                     bfe7+
                     age_group +income +edu+
                     survey_2006+
                     survey_2010+
                     survey_2016,
                   data=subset(data),
                   na.action = na.exclude)
summary(bribery_documents)
bribery_documents_cov <- vcovHC(bribery_documents, type = "HC")
bribery_documents_se <- sqrt(diag(bribery_documents_cov))


bribery_documents_b1945<-lm(bribery_documents~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=data_before_1945,
                            na.action = na.exclude)
summary(bribery_documents_b1945)
bribery_documents_b1945_cov <- vcovHC(bribery_documents_b1945, type = "HC")
bribery_documents_b1945_se <- sqrt(diag(bribery_documents_b1945_cov))


bribery_documents_1945_1990<-lm(bribery_documents~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=data_before_1945_1990,
                            na.action = na.exclude)
summary(bribery_documents_1945_1990)
bribery_documents_1945_1990_cov <- vcovHC(bribery_documents_1945_1990, type = "HC")
bribery_documents_1945_1990_se <- sqrt(diag(bribery_documents_1945_1990_cov))


bribery_documents_1990_2016<-lm(bribery_documents~treat+
                                  POINT_X+POINT_Y+zagreb_distance+
                                  bfe1+bfe2+bfe3+
                                  bfe4+bfe5+bfe6+
                                  bfe7+
                                  age_group +income +edu+
                                  survey_2006+
                                  survey_2010+
                                  survey_2016,
                                data=data_before_1990_2016,
                                na.action = na.exclude)
summary(bribery_documents_1990_2016)
bribery_documents_1990_2016_cov <- vcovHC(bribery_documents_1990_2016, type = "HC")
bribery_documents_1990_2016_se <- sqrt(diag(bribery_documents_1990_2016_cov))



stargazer(bribery_documents_b1945, 
          bribery_documents_1945_1990,
          bribery_documents_1990_2016,
          keep = c("treat"),
          type = "text",
          covariate.labels=c("Military Territory"),
          se=list(bribery_documents_b1945_se, 
                  bribery_documents_1945_1990_se, 
                  bribery_documents_1990_2016_se),
          dep.var.labels.include = F,
          omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE)

################
#bribery_education#
################
bribery_social_sec<-lm(bribery_social_sec~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=subset(data),
                       na.action = na.exclude)
summary(bribery_social_sec)
bribery_social_sec_cov <- vcovHC(bribery_social_sec, type = "HC")
bribery_social_sec_se <- sqrt(diag(bribery_social_sec_cov))


bribery_social_sec_b1945<-lm(bribery_social_sec~treat+
                               POINT_X+POINT_Y+zagreb_distance+
                               bfe1+bfe2+bfe3+
                               bfe4+bfe5+bfe6+
                               bfe7+
                               age_group +income +edu+
                               survey_2006+
                               survey_2010+
                               survey_2016,
                             data=data_before_1945,
                             na.action = na.exclude)
summary(bribery_social_sec_b1945)
bribery_social_sec_b1945_cov <- vcovHC(bribery_social_sec_b1945, type = "HC")
bribery_social_sec_b1945_se <- sqrt(diag(bribery_social_sec_b1945_cov))


bribery_social_sec_1945_1990<-lm(bribery_social_sec~treat+
                                   POINT_X+POINT_Y+zagreb_distance+
                                   bfe1+bfe2+bfe3+
                                   bfe4+bfe5+bfe6+
                                   bfe7+
                                   age_group +income +edu+
                                   survey_2006+
                                   survey_2010+
                                   survey_2016,
                                 data=data_before_1945_1990,
                                 na.action = na.exclude)
summary(bribery_social_sec_1945_1990)
bribery_social_sec_1945_1990_cov <- vcovHC(bribery_social_sec_1945_1990, type = "HC")
bribery_social_sec_1945_1990_se <- sqrt(diag(bribery_social_sec_1945_1990_cov))



bribery_social_sec_1990_2016<-lm(bribery_social_sec~treat+
                                   POINT_X+POINT_Y+zagreb_distance+
                                   bfe1+bfe2+bfe3+
                                   bfe4+bfe5+bfe6+
                                   bfe7+
                                   age_group +income +edu+
                                   survey_2006+
                                   survey_2010+
                                   survey_2016,
                                 data=data_before_1990_2016,
                                 na.action = na.exclude)
summary(bribery_social_sec_1990_2016)
bribery_social_sec_1990_2016_cov <- vcovHC(bribery_social_sec_1990_2016, type = "HC")
bribery_social_sec_1990_2016_se <- sqrt(diag(bribery_social_sec_1990_2016_cov))


###############
#part_demonstr#
###############

part_demonstr_20km<-lm(part_demonstr~treat+
                    POINT_X+POINT_Y+zagreb_distance+
                    bfe1+bfe2+bfe3+
                    bfe4+bfe5+bfe6+
                    bfe7+
                    age_group +income +edu+
                    survey_2006+
                    survey_2010+
                    survey_2016,
                  data=subset(data, krajna6_distance<20), na.action=na.exclude)
summary(part_demonstr_20km)


part_demonstr_40km<-lm(part_demonstr~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=subset(data, krajna6_distance<40), na.action=na.exclude)
summary(part_demonstr_40km)


part_demonstr_60km<-lm(part_demonstr~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=subset(data, krajna6_distance<60), na.action=na.exclude)
summary(part_demonstr_60km)

part_demonstr<-lm(part_demonstr~treat+
                     POINT_X+POINT_Y+zagreb_distance+
                     bfe1+bfe2+bfe3+
                     bfe4+bfe5+bfe6+
                     bfe7+
                     age_group +income +edu+
                     survey_2006+
                     survey_2010+
                     survey_2016,
                    data=subset(data), na.action=na.exclude)
summary(part_demonstr)
part_demonstr_cov <- vcovHC(part_demonstr, type = "HC")
part_demonstr_se <- sqrt(diag(part_demonstr_cov))

part_demonstr_b1945<-lm(part_demonstr~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=data_before_1945, na.action=na.exclude)
summary(part_demonstr_b1945)
part_demonstr_b1945_cov <- vcovHC(part_demonstr_b1945, type = "HC")
part_demonstr_b1945_se <- sqrt(diag(part_demonstr_b1945_cov))


part_demonstr_1945_1990<-lm(part_demonstr~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=data_before_1945_1990, na.action=na.exclude)
summary(part_demonstr_1945_1990)
part_demonstr_1945_1990_cov <- vcovHC(part_demonstr_1945_1990, type = "HC")
part_demonstr_1945_1990_se <- sqrt(diag(part_demonstr_1945_1990_cov))

part_demonstr_1990_2016<-lm(part_demonstr~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=data_before_1990_2016, na.action=na.exclude)
summary(part_demonstr_1990_2016)
part_demonstr_1990_2016_cov <- vcovHC(part_demonstr_1990_2016, type = "HC")
part_demonstr_1990_2016_se <- sqrt(diag(part_demonstr_1990_2016_cov))


###########
#sign_petition#
###############

sign_petition_20km<-lm(sign_petition~treat+
                    POINT_X+POINT_Y+zagreb_distance+
                    bfe1+bfe2+bfe3+
                    bfe4+bfe5+bfe6+
                    bfe7+
                    age_group +income +edu+
                    survey_2006+
                    survey_2010+
                    survey_2016,
                  data=subset(data, krajna6_distance<20), na.action = na.exclude)
summary(sign_petition_20km)


sign_petition_40km<-lm(sign_petition~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=subset(data, krajna6_distance<40), na.action = na.exclude)
summary(sign_petition_40km)


sign_petition_60km<-lm(sign_petition~treat+
                         POINT_X+POINT_Y+zagreb_distance+
                         bfe1+bfe2+bfe3+
                         bfe4+bfe5+bfe6+
                         bfe7+
                         age_group +income +edu+
                         survey_2006+
                         survey_2010+
                         survey_2016,
                       data=subset(data, krajna6_distance<60), na.action = na.exclude)
summary(sign_petition_60km)

sign_petition<-lm(sign_petition~treat+
                     POINT_X+POINT_Y+zagreb_distance+
                     bfe1+bfe2+bfe3+
                     bfe4+bfe5+bfe6+
                     bfe7+
                     age_group +income +edu+
                     survey_2006+
                     survey_2010+
                     survey_2016,
                data=subset(data), na.action = na.exclude)
summary(sign_petition)
sign_petition_cov <- vcovHC(sign_petition, type = "HC")
sign_petition_se <- sqrt(diag(sign_petition_cov))


sign_petition_b1945<-lm(sign_petition~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=data_before_1945, na.action = na.exclude)
summary(sign_petition_b1945)
sign_petition_b1945_cov <- vcovHC(sign_petition_b1945, type = "HC")
sign_petition_b1945_se <- sqrt(diag(sign_petition_b1945_cov))


sign_petition_1945_1990<-lm(sign_petition~treat+
                          POINT_X+POINT_Y+zagreb_distance+
                          bfe1+bfe2+bfe3+
                          bfe4+bfe5+bfe6+
                          bfe7+
                          age_group +income +edu+
                          survey_2006+
                          survey_2010+
                          survey_2016,
                        data=data_before_1945_1990, na.action = na.exclude)
summary(sign_petition_1945_1990)
sign_petition_1945_1990_cov <- vcovHC(sign_petition_1945_1990, type = "HC")
sign_petition_1945_1990_se <- sqrt(diag(sign_petition_1945_1990_cov))


sign_petition_1990_2016<-lm(sign_petition~treat+
                              POINT_X+POINT_Y+zagreb_distance+
                              bfe1+bfe2+bfe3+
                              bfe4+bfe5+bfe6+
                              bfe7+
                              age_group +income +edu+
                              survey_2006+
                              survey_2010+
                              survey_2016,
                            data=data_before_1990_2016, na.action = na.exclude)
summary(sign_petition_1990_2016)
sign_petition_1990_2016_cov <- vcovHC(sign_petition_1990_2016, type = "HC")
sign_petition_1990_2016_se <- sqrt(diag(sign_petition_1990_2016_cov))



#############
#Final Table#
#############
list_models<-list(trust_people,
                  trust_family,
                  trust_foreign_investors,
                  trust_presidency,
                  bribery_road_police,
                  bribery_unempl,
                  part_demonstr,
                  sign_petition)

public_opinion<-stargazer(list_models,
          object.names = F,
          column.labels= c("\\shortstack{Trust People,\\\\ 2006-16}", 
                           "\\shortstack{Trust Family,\\\\ 2016}", 
                           "\\shortstack{Trust Foreign\\\\ Investors,\\\\ 2006-16}",
                           "\\shortstack{Trust\\\\ Presidency,\\\\ 2006-16}",
                           "\\shortstack{Bribery\\\\ Road Police,\\\\ 2006-16}",
                           "\\shortstack{Bribery \\\\ Unemployment,\\\\ 2006-16}",
                           "\\shortstack{Participation\\\\ Demonstrations,\\\\ 2006-16}",
                           "\\shortstack{Sign Petition,\\\\ 2016-16}"),
          keep = c("treat"),
          covariate.labels=c("Military Territory"),
          model.names = F,
          se=list(trust_people_se,
                  trust_family_se,
                  trust_foreign_investors_se,
                  trust_presidency_se,
                  bribery_road_police_se,
                  bribery_unempl_se,
                  part_demonstr_se,
                  sign_petition_se),
          dep.var.caption = "Dependent variable:",
          dep.var.labels.include = F,
          omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
          add.lines=list(c("Mean", round(mean(data$trust_people, na.rm=T), digits = 3), 
                           round(mean(subset(data, survey_2016==1)$trust_family, na.rm=T), 3),
                           round(mean(data$trust_foreign_investors, na.rm=T), 3),
                           round(mean(data$trust_presidency, na.rm=T), 3),
                           round(mean(data$bribery_road_police, na.rm=T), 3),
                           round(mean(data$bribery_unempl, na.rm=T), 3),
                           round(mean(data$part_demonstr, na.rm=T), 3),
                           round(mean(data$sign_petition, na.rm=T), 3)),
                         c("SD", round(sd(data$trust_people, na.rm=T), digits = 3), 
                           round(sd(subset(data, survey_2016==1)$trust_family, na.rm=T), 3),
                           round(sd(data$trust_foreign_investors, na.rm=T), 3),
                           round(sd(data$trust_presidency, na.rm=T), 3),
                           round(sd(data$bribery_road_police, na.rm=T), 3),
                           round(sd(data$bribery_unempl, na.rm=T), 3),
                           round(sd(data$part_demonstr, na.rm=T), 3),
                           round(sd(data$sign_petition, na.rm=T), 3)),
                         c("Boundary FE", rep("Yes", 9)),
                         c("Survey Year FE", rep("Yes", 9)),
                         c("Demographic Covariates", rep("Yes", 9))))
          
note_latex <- "\\multicolumn{9}{l} {\\parbox[t]{29cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is the individual. Data is based on EBRD LITS 2006, 2010 and 2016.
The trust variables range 
from 1 (complete distrust) to 5 (complete trust). The bribery variables range from 1 (never) to 5 (always). 
Participation in demonstration and signing a 
petition ranges from 0 (never and might) to 1 (have done). Every regression includes the following covariates:
age, perception of income, and education.
See appendix for details on variables and data sources.}}"
public_opinion [grepl("Note", public_opinion)] <- note_latex
cat(public_opinion, file="./Dropbox/Legacies_Project/Paper/tables/table_4.tex", sep="\n")


return_graphs_coef_pres<-function(models_run, models_run_alias){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  irislabs1 <- rev(big_data$dv_official)
  irislabs2 <- rev(big_data$no_obs)
  #return(big_data)
  mean_feature <- ggplot(data = big_data, aes(x=estimate, y=as.numeric(reorder(dv_official, desc(counter))))) +
    geom_point(size = 2) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs1),
                       labels = irislabs1,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    #scale_y_discrete(name = "") +
    scale_x_continuous(name = "Estimate") +
    #ggtitle(glue("{title_spec}")) +
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=18),
          axis.text.y = element_text(size=18),
          axis.title=element_text(size=18))  
  #return(big_data)
  return(mean_feature)
}

#By distances
return_graphs_coef_band<-function(models_run, models_run_alias, y_label){
  #Process:
  #-takes list of regression objects, associates them with a 
  #list of "official" label, keep the "treat" variable and
  #produce ggplots with those
  
  #Input:
  #-List of Regression model object after lm
  #-List of label for the DV that you want displayed in ggplot
  
  #Output:
  #-A ggplot with regression coefficients and robust CI
  
  names(models_run)<-models_run_alias
  df<-list()
  length(df)<-length(models_run)
  counter = 1
  for (i in 1:length(models_run)){
    #print(i)
    x<-models_run[[i]]
    name_x<-names(models_run)[i]
    #print(name_x)
    model_tidy<-rerun_regression(x, name_x)
    model_tidy$counter<-counter
    #print(model_tidy)
    df[[i]] = model_tidy
    counter = counter +1}
  print(df)
  big_data = do.call(rbind, df)
  irislabs1 <- big_data$dv_official
  irislabs2 <- big_data$no_obs
  #return(big_data)
  mean_feature <- ggplot(data = big_data, aes(x=estimate, y=as.numeric(reorder(dv_official, counter)))) +
    geom_point(size = 2) +
    geom_errorbar(aes(xmin = lb_95, xmax = ub_95), size = 0.9, width = 0)+
    scale_y_continuous(name = "",
                       breaks = 1:length(irislabs1),
                       limits = c(0.5,length(irislabs1)+0.5),
                       labels = irislabs1,
                       sec.axis = sec_axis(~.,
                                           breaks = 1:length(irislabs2),
                                           labels = irislabs2))+
    geom_vline(xintercept = 0) +
    theme_bw() +
    #scale_y_discrete(name = "") +
    scale_x_continuous(name = glue("{y_label}")) +
    #ggtitle(glue("{title_spec}")) +
    coord_flip()+
    theme(plot.title = element_text(hjust = 0.5),
          axis.text.x = element_text(size=18),
          axis.text.y = element_text(size=18),
          axis.title=element_text(size=18))  
  #return(big_data)
  return(mean_feature)
}



models_run<-list(trust_people_20km,
                 trust_people_40km,
                 trust_people_60km)


column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

library("glue")
sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Others")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18a.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Trust Generations
models_run<-list(trust_people_b1945,
                 trust_people_1945_1990,
                 trust_people_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Others")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19a.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Trust Family
models_run<-list(trust_family_20km,
                 trust_family_40km,
                 trust_family_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Family")

ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18b.jpg', 
       height=10, width=18, units = "cm", dpi=300)


models_run<-list(trust_family_b1945,
                 trust_family_1945_1990,
                 trust_family_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)


sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Family")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19b.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Trust Foreign Investors

models_run<-list(trust_foreign_investors_20km,
                 trust_foreign_investors_40km,
                 trust_foreign_investors_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Foreign Investors")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18c.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Trust foreign Investors Generations
models_run<-list(trust_foreign_investors_b1945,
                 trust_foreign_investors_1945_1990,
                 trust_foreign_investors_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Foreign Investors")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19c.jpg', 
       height=10, width=18, units = "cm", dpi=300)

#Trust Presidency
models_run<-list(trust_presidency_20km,
                 trust_presidency_40km,
                 trust_presidency_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Presidency")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18d.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Trust Presidency
models_run<-list(trust_presidency_b1945,
                 trust_presidency_1945_1990,
                 trust_presidency_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Trust Presidency")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19d.jpg', 
       height=10, width=18, units = "cm", dpi=300)

#Bribery
models_run<-list(bribery_road_police_20km,
                 bribery_road_police_40km,
                 bribery_road_police_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)


sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Bribery Road Police")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18e.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Bribery Generation
models_run<-list(bribery_road_police_b1945,
                 bribery_road_police_1945_1990,
                 bribery_road_police_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Bribery Road Police")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19e.jpg', 
       height=10, width=18, units = "cm", dpi=300)


models_run<-list(bribery_unempl_20km,
                 bribery_unempl_40km,
                 bribery_unempl_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Bribery Unemployment")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18f.jpg', 
       height=10, width=18, units = "cm", dpi=300)

#Unemployment Generations
models_run<-list(bribery_unempl_b1945,
                 bribery_unempl_1945_1990,
                 bribery_unempl_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Bribery Unemployment")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19f.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Participate Demonstrations
models_run<-list(part_demonstr_20km,
                 part_demonstr_40km,
                 part_demonstr_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Participation in Demonstrations")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18g.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Participate Demonstrations
models_run<-list(part_demonstr_b1945,
                 part_demonstr_1945_1990,
                 part_demonstr_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)

sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Participation in Demonstrations")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19g.jpg', 
       height=10, width=18, units = "cm", dpi=300)


#Sign petitions
models_run<-list(sign_petition_20km,
                 sign_petition_40km,
                 sign_petition_60km)

column_labels= c("\\shortstack{20km}", 
                 "\\shortstack{40km}",
                 "\\shortstack{60km}")

models_run_alias_clean<-clean_col_labels(column_labels)


sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Sign Petitions")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A18h.jpg', 
       height=10, width=18, units = "cm", dpi=300)

#Sign petitions
models_run<-list(sign_petition_b1945,
                 sign_petition_1945_1990,
                 sign_petition_1990_2016)

column_labels= c("\\shortstack{Born\\\\\ Before\\\\\ 1945}", 
                 "\\shortstack{Born\\\\\ Between\\\\\ 1945-1990}",
                 "\\shortstack{Born\\\\\ Between\\\\\ 1990-2016}")

models_run_alias_clean<-clean_col_labels(column_labels)


sdsd<-return_graphs_coef_band(models_run, models_run_alias_clean, "Sign Petitions")
ggsave(sdsd, file = './Dropbox/Legacies_Project/Paper/figures/Figure_A19h.jpg', 
       height=10, width=18, units = "cm", dpi=300)




#####################
#Effect of Communism#
#####################

#Torture against oneself
q925a<-lm(q925a~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=subset(data, !is.na(data$q925a)))
summary(q925a)

#Torture against family
q925c<-lm(q925c~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=subset(data, !is.na(data$q925c)))
summary(q925c)

#Torture against relatives
q925d<-lm(q925d~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=data)
summary(q925d)


##############
#Final Table#
#############


relevant_vars<-c("q925a",
                 "q925d",
                 "q925c")


list_models<-list(q925a,
                  q925d,
                  q925c)

column_labels<-c("\\shortstack{Torture\\\\ against Oneself}", 
                 "\\shortstack{Torture\\\\ against Family}",
                 "\\shortstack{Torture against\\\\ Grandparents}")

means<-list()
sds<-list()
for (i in relevant_vars){
  item_mean<-round(mean(data[[i]], na.rm=T), 3)
  item_sd<-round(sd(data[[i]], na.rm=T), 3)
  means[[length(means) + 1]] <- item_mean 
  sds[[length(sds) + 1]] <- item_sd 
}
means<-unlist(means)
sds<-unlist(sds)



eff_com=stargazer(list_models,
                  object.names = F,
                  column.labels= column_labels,
                  keep = c("treat"),
                  covariate.labels=c("Military Territory"),
                  dep.var.caption = "Dependent variable:",
                  dep.var.labels.include = F,
                  omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                  add.lines=list(c("Mean", means),
                                 c("SD", sds),
                                 c("Boundary FE", rep("Yes", 9)),
                                 c("Survey Year FE", rep("Yes", 9)),
                                 c("Demographic Covariates", rep("Yes", 9))))
note_latex <- "\\multicolumn{4}{l} {\\parbox[t]{15cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from Logit regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is the individual. 
See codebook for data sources and how the variables were calculated.}}"
eff_com [grepl("Note", eff_com)] <- note_latex
cat(eff_com, file="./Dropbox/Legacies_Project/Paper/tables/table_A13.tex", sep="\n")


#############
#Effect WW2#
#############

#Were you, your parents or any of your grandparents physically injured or 
#were your parents or any of your grandparents killed during the Second World War?
#0  Not mentioned;  1  Mentioned"

q924a<-lm(q924a~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=subset(data, data$q924a!=99))
summary(q924a)

#Did you, your parents or any of your grandparents have to move as a result of the Second World 
q924b<-lm(q924b~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=subset(data2))
summary(q924b)


##########################
#Effect War in Yugoslavia#
##########################

#Were you or any member of your household physically injured as a result 
#of the conflict in Yugoslavia (1991-2001) "
#0  Not mentioned;  1  Mentioned"
q924c<-lm(q924c~treat+
            POINT_X+POINT_Y+zagreb_distance+
            bfe1+bfe2+bfe3+
            bfe4+bfe5+bfe6+
            bfe7+
            age_group +income +edu+
            survey_2006+
            survey_2010+
            survey_2016,
          data=subset(data))
summary(q924c)

#Did your household have to move as a result
#of the conflict in Yugoslavia (1991-2001) 

relevant_vars<-c("q924a", "q924b", "q924c")
list_models<-list(q924a, q924b, q924c)

column_labels= c("\\shortstack{Injury\\\\ during WW2}", 
                 "\\shortstack{Had to Move\\\\ during WW2}",
                 "\\shortstack{Injury during\\\\ Yugoslav Conflict}")

means<-list()
sds<-list()
for (i in relevant_vars){
  item_mean<-round(mean(data[[i]], na.rm=T), 3)
  item_sd<-round(sd(data[[i]], na.rm=T), 3)
  means[[length(means) + 1]] <- item_mean 
  sds[[length(sds) + 1]] <- item_sd 
}
means<-unlist(means)
sds<-unlist(sds)


eff_wars =stargazer(list_models,
                    object.names = F,
                    column.labels= column_labels,
                    keep = c("treat"),
                    covariate.labels=c("Military Territory"),
                    dep.var.caption = "Dependent variable:",
                    dep.var.labels.include = F,
                    omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                    add.lines=list(c("Mean", means),
                                   c("SD", sds),
                                   c("Boundary FE", rep("Yes", 9)),
                                   c("Survey Year FE", rep("Yes", 9)),
                                   c("Demographic Covariates", rep("Yes", 9))))
note_latex <- "\\multicolumn{4}{l} {\\parbox[t]{15cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from Logit. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. The unit of analysis is the individual. 
See codebook for variable descriptions.}}"
eff_wars [grepl("Note", eff_wars)] <- note_latex
cat(eff_wars, file="./Dropbox/Legacies_Project/Paper/tables/table_A15.tex", sep="\n")

